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DECUS Program Library Write-up 


DECUS NO. 8-590 


MATRIX INVERSION 


ABSTRACT. 

* 

This subroutine inverts a matrix of reai numbers residing in core 
in floating format (three words). The Gauss-Jordan method of elimination is used 
to perform the inversion. 


REQUIREMENTS 

Storage . 

This subroutine occupies 737g + 2.N0 locations anywhere in core 
and four locations in page 0. N0 is the matrix order. The beginning of the sub¬ 
routine must be at a page boundary. The listing which accompanies this write-up 
shows the subprogram beginning at location 200 8 - This subroutine uses the basic 
Floating Point Package (F.P.P.) and a patch of 7 instructions (RFSGE) which 
correspond to an additional operation in F.P.P. and must be in the same field 
that the F.P.P. 

If F.P.P. is not in the same memory field than the routine one, 
the extended memory patch (DEC 8 364 ) must be used. 

Equipment : PDP8 I or E with EAE or 

PDP8 L with instruction set simulator (DEC 8 - 17U). 

USAGE OF THE SUBROUTINE. 

Loading. 

The PAL symbolic tape concerns the subroutine only, and may be 
assembled into an alternate location. (Warning : don't forget to kee p free 2N0 
locations after the end of the subroutine). 










Assembling parameter s. 


Before the assembling, let set up the following parameters : 

FPPFLD = X where X is equal to 00, 10, 20, _ etc according 

to the F.P.P. resident field 0, 1 , 2 etc.... 


CURFLD = Y where Y is equal to 00, 10, 20,_ etc according 

to the subroutine, the calling program and data memory field (the same for 
the 3). 


Calling sequence . 



The elements a^ of the matrix must be in order a , a 
.. a nn * (i.e. row after row). 



Initialize : 

ADR0,... / first address of the matrix array 

N0,.. / matrix order 

JMS (INVMAT) 

Return will be to PC+l with the inverse of matrix 
in place of the original matrix array and the determinant is returned as the 
value of DET. 


Utilization of internal subroutine 


s. 



Some internal subroutines may be used for deposing data into the 
matrix array, or picking up results from the same array. 

Warning : The first element is a and the last a 

oo n-1’ n-l 

(instead of a^ to a^ ^ usually used). 

- Deposing data : Initialize II and J 

JMS (RCALIJ.) 

Return will be to PC+l with the address of the element (a. .) in 
PTAR (the address of the exponent). 
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In the same way, we can use for an element a, . the routine 

• • , j 

RCAL KJ after initializing K and J. The address of the element will be in PTAIJ 

A third subroutine RCALIK is disponible for an element a. , . The 
,, . i,k 

address is m PTAIJ. 


“ Picking up data : As for deposing data. 


- Utilization of loops tests : If INVMAT has not still 
been used, initialize N = N0 - 1. 


Test : J = N ? 

PC JMS R TEST j / J = N ? 

PC + 1 RETURN FROM TEST / j / N 

PC + 2 RETURN FROM TEST / J = N 

Test : II = N ? 


PC JMS R TEST I / II = N ? 

PC + 1 RETURN FROM TEST / II f N 

PC + 2 RETURN FROM TEST / II f N 


RESTRICTIONS. 


field. 


The calling program, the subroutine and the data must be in the same 


If we use an other version than YQYA-PB for the F.P.P., let's 
verify that the entry point (for additional pseudo-operator FSGE - ll in F.P.P.) 
of RFSGE may be in address 6555 of TABLE 6 of the F.P.P. One can change the code 
of FSGE and the address entry point in TABLE 6 of the F.P.P. 


USAGE OF THE PROGRAM (BINARY PAPER TAPE). 


Loading . 

A routine which allows input of data and output of results on T.T.Y. 

is supplied as a binary tape which may be loaded with the BIN loader after the 
F.P.P. 
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- Operating procedure . 


. Starting address : 1000 
. Respond to question 


MATRIX ORDER : 

. Type values of coefficients row after row. Each value will 
be separated by a space. 

. Respond to question. 

CHANGE VALUE ? (0 = No : 1 - yes) ^ 

. The coefficient of the inverse matrix are printed out row 

after row. 

. Respond to question : INVERSE INVERSE ? (0 = No • l = yes ) 


- Informations. 


. Time of execution (see figure). 

. Example of utilization with a matrix of which value 
determinant is near 0 . 


of 



We use the pseudo F.P.P. operator- INPUT = 13 and OUTPUT = 14, 
and the F.P.P. used must be accomoded for (in TABLED). 
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u 

3 H 


2 - 


1 - 



EXECUTION TIME 
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MATRIX ORDER:4 


3 8 5 

3 8 5 

4 6 8 

4 2 9 

CHANGE VALUE70 
INVERSE MATRIX 

+0.3020882E+04 
- 0 »402250 7 E+0 3 
-0.1408877E+04 
+0.2001247E+03 
INVERSE INVERSE?1 
+0.2998024E+01 
+0.2998026E+01 
+0.3996766E+01 
+0.3996387E+01 
MATRIX ORDER:4 
3 5 2 

6 8 3 

45 7 9 

5 3 9 

CHANGE VALUE? 1 
ROW: 3 
COLUMN: 1 
VALUE: 8 

CHANGE VALUE?0 
INVERSE MATRIX 

-0.5833323E+00 
+0.4999994E+00 
+0.I666664E+00 
-0.8333335E-01 
INVERSE INVERSE?0 


6 

5.995 

8 

7 

-0 » 3001876E+04 
+ 0.4002499 E+0 3 
+0.14008 7 4E + 04 
-0*2001242E+03 

+ 0.799461 1E + 0 1 
+0.7994611E+01 
+0.5991228E+01 
+0.1990233E+01 

1 

12 

4 

1 


-0.2247471E+00 
+0.7575742E-01 
+ 0.858 58 50E-01 
+ 0.1237 37 3E + 00 


-0.3101118E+02 
+0.3501489E+01 
+ 0.1300 521E+02 
-0.7400512E-03 

+0.4996701E+01 
+ 0•499 6742E + 01 
+0•7994650E+01 
+0.8994015E+01 


+0.1025250E+01 
-0.4242416E+00 
- 10 . 41 41 40 9E + 00 
-0*1262625E+00 


+0*1700648E+02 
-0.2000865E+01 
-0.7003027E+01 
+0.4293918E-03 

+0.5996038E+01 
+0.3991041E+01 
+ 0.79 9 3520E + 01 
+0.6992763E+01 


-0.8207057t+00 
+0.28787818+00 
+0.4595953E+00 
+0.1035352E+00 


MATRIX ORDER 
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R PIP 


TTY:<INVMAT 


/ROUTINE D'INVERSION DE MATR1CE 
FPPFLD=00 /RESIDENT FIELD OF iHE F.P.P. 

CURFLD=00 /FIELD WHERE ARE THE MATRIX & THE NAT• INV. ROUTINE 
/ADR0 * CONTAIN THE 1ST AUKESSE OF TriE MATRIX i 
/N0 , THE ORDER OF THE MATRIX 
/ 

/THE ROUTINE USES A BUFFER EQUAL TO 2*N0 JUST AF1EK IIS LAST 
/ ADRESSE OR AFTER THE I/O SUBPROGRAMME LAST ONE #IF IN CORE 
/CALLING SEQUENCE 
/JMS I (IwVrtAT) 

/ RETURN FROM THE ROUTINE WITH THE INVERTED MATRIX 
/ AT THE SAME PLACE OR WITH DE1=0 1 F 

/ THE DETERMINANT IS NUL 

/DATA FIELD IS RESET TO (CALFLD) i AC=0 
/ 

FIELD 0 /SET TO 0#1>2«. ACCORDING TO FPP RESIDENT FIELD 
*7 

S600 

*6bSb 

RFSGE 

/ 

/ 

*4400 /EVERY FREE PLACE IS GOOD 
RFSGE> HLT /REAL I EE A SKIP IF FAC > OR = TO 0 
CL A 

TAD 4b 
SMA CLA 
I SZ I Ab6bb 
JMP I RFSGE 
Ab6bb»b6bb 
/ 

FSGE=11 
/ 

FIELD 0 / SET TO 0«1«2».. ACCORDING TO CCUKFLD> 

*70 /CONSTANTS IN EERO PAGE NOT NECESSARLY AT C70> 

PTAI J,0 
PTAR*0 
PTJK*0 
PTIK,0 

*200 f 

INVMAT* HL T 
CIF+FPPFLD 
JMS I 7 
F GET FUN 
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FPUT DET /DETERMINANT^1. 

FLAT 

TAD PTIK0 
TAD N0 
DCA PTJK0 
CL A Cm A 
TAD i>i0 
DCA N 
JMS INVER 

JMP DETNUL /RLTOUR SI DETERMiNANT=0 
DCA L /RESTAURLR L’ORDRL DL LA NAT RICE 
RESTAU* TAD L 

CMA IAC / K=N-L+l 

TAD N / UN DOIT METTRE N-L CAR ON PART A L = 0 

DCA K 

TAD PTIK0 

TAD K 

DCA PT IK 

TAD I FTIK 

DCA J / J= I K< K) 

JMS RJMK /J-K 
SPA SNA CLA 
JMP .+7 
DCA II /l=0 
JMS RCALIK 
JMS RCAL1J 

JMS RSAVAA /HOLD=A(i#K) l A(I#K)=A(I#j) ; A<I » J) = -rtULU 
JMS RTESTI /I=N ? 

JMP .-4 /NON 1 = 1 + 1 I* AIT DANS TESII 

TAD PTJK0 

TAD K 

DCA PTJK 

TAD I PTJK 

DCA II /1=JKC K > 

JMS RIMK /1 - K 
SPA SNA CLA 
JMP .+7 
DCA J //J=0 
JMS RCALKJ 
JMS RCAlIJ 

JMS RSAVAA /HOLD=A(K>J) > AC K* J > =-A(1 * J) JA(1>J>=SAVE 
JMS RTESTJ /J=N ? 

JMP .-4 /NON 

TAD K /Si K = 0 J L=N 

SNA CLA /L=N ? 

JMP I INVMAT /OU1 EXIT FIN DL U’iNVLRSION 
IS<i L 

JMP RLSTAU /NON 

/ 

DETNUL> DCA DET 
DCA DET+1 
DCA DLT + 2 
JMP 1 INVMAT 

RIMK* HLT /CACUL DL 1-K 
TAD K 
CIA 

TAD il 

JMP I RIMK /RESULT DANS AC 
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RJMK# HLT /CALCUL DE J-K 
TAD K 
CIA 
TAD J 

JMP I RJMK 


RCALIJ^HLT /CALCUL DE L‘ADRESSE DE ( A(i*J> > 

TAD 11 

JMS RMULTN / I *N 
TAD J 

JMS MUY3P / (I (I*N)+J>*3 ♦ADRESSE DEPART 
DCA PTAR 
JMP I RCALIJ 

/ 

RCALIKjHLT /CALCUL DE L‘ADRESSE DE < AC I *K) ) 

TAD II 
JMS RMULTN 
TAD * 

JMS MUY3P 
DCA PTAIJ 
JMP I RCAL1K 

/ 

RCALKJ^HLT / POUR ELEMENT C A(K,J> ) 

TAD K 

JMS RMULTN 
TAD J 
JMS MUY3P 
DCA PTAIJ 
JMP I RCALKJ 

/ 

MUY3P*HLT/ NB PRIS DANS AC *3 +ADRESSE DEPART(ADR0> 

MQL MUY 
3 

MQA 

TAD ADRfe) 

JMP I MUY3P /RETOUR AVEC RESULTAT DANS AC 

/ 

RTESTN* HLT 

TAD I RTESTN /RECUPEREK L‘ARGUMENT SUiVANT LE JMS 
CM A I AC 
TAD N 

SNA CLA /ARG=N ? 

JMP •+3 /OUI 

1SZ I RTESTN/INCREMENTER L‘ARGUMENT 
ISZ RTESTN /♦1ADRESSE RETUUR 
ISZ RTESTN / 

JMP I RTESTN/RETOUR A ARG+1 SI ARG=NJ ARG+2 Si ARG DIP- DE N 


9 



/ 

/ 

/ 

RMUL TN# HLT 
MQL MUY 
N0> 3 
MQA 

JMP I RMULTN 

/ 

ADR0* 1000 
N» 2 
L* 0 

PTIA0* DEBIA 
PTJA0*0 
DET* 0 J 0 J 0 
/ 

PAGE 

/ 

/ 

/routine o'inversion ppt due 

INVER* HLT 
DCA K 

A1 * DCA AMAX+1 
DCA AMAX+2 

TAD A /RECHERCHE DE L 'ELENENT MAX UE CHAwUE LlGNE ONREPERE 
DCA II / LIGNE ET COL DE L * ELEN ENT MAX D A»nS 1ACA) t.T JACK) 
LI GSUI * TAD A 
DCA J /J=K 
COLSIJI * JMS RCALIJ 
TAD PTIK0 
TAD K 

DCA PTIK /POINTEUR SUR IACK) 

TAD PTJK0 
TAD K 

DCA PTJA / ” •* JACK) 

JMS ABSVL / CAlCUL VAL. ABS. DE AMAX ET AC1#J) 

JMS KTESTJ /J=N ? 

JMP COLSUI /NON 
JMS RTESTI /I -N ? 

JMP LIGSUI /NON 

CIF+FPPFLD 

JMS I 7 

FGET AMAX 

FEXT 

TAD 4S 

SNA CLA 

JMP I INVER /OU1 EXIT 

TAD 1 PTIK 

DCA II /I=IACA) 

JMS RiMK /l-A 
SPA SNA CLA 
JMP JJKK /I-A =0 
DCA J /1~K >0 
JMS RCALAJ 
JMS RCALIJ 

JMS RHOLAA /HULD=ACA*J> l A C A* J ) = A C 1 , J ) i AC 1 * J ) = -HULL) 

JMS KTESTJ /J=N ? 

JMP .-A 
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JJKK# TAD I PTJK 
DCA J /J= JK(K) 

JMS RJMK 

SPA SNA CLA 

JMP ASSURM 

DCA II 

JMS RCALIK 

JMS RCALIJ 

JMS RHOLAA 

JMS RTEST1 /1 = t\l ? 

JMP .-4 

ASSURM#DCA 11 /i=0 
JMS RIMK 

SNA CLA /1-K = 0 ? 

JMP .+11 
JMS RCALIK 
CIF+FPPFLD 
JMS I 7 
FGET I PTAIJ 
FMPY FM1 
FDIV AMAX 

FPUT I PTAIJ / A (I # K ) =-A(l»K)/ArtAX 
FEXT 

JMS RTEST1 /1 = N ? 

JMP ASSURM+1 
APASA# DCA II 
DCA J 

JMS RIMK /ATTENTION IL Y A 1 JMP APASA+1 Ef + 2 
SNA CLA 
JMP FINBOU 
JMS RJMK 

SZA CLA /J-K = 0 ? 

JMS APAMA /NON 
FINBOU#JMS RTESTJ 
JMP APASA+2 
JMS RTESTI 
JMP APASA+1 
DCA J 

DI VAMA# JMS RJMK 
SNA CLA 
JMP .+10 
JMS RCALKJ 
CIF+FPPFLD 
JMS I 7 
FGET I PTAIJ 
FDIV AMAX 
FPUT I PTAIJ 
FEXT 

JMS RTESTJ 
JMP DIVAMA 
TAD K 

JMS RMULTN 
TAD K 
JMS MUY3P 
DCA PTAR 
JMS DETAMA 
JMS RTESTN 
K# 0 
SKP 
JMP A 1 
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I SZ INVER 
JMP I INVER 


RTESTJjHLT 
JMS RTESTN 
J>0 

ISZ RTESTJ 
JMP I RTESTJ 

/ 

RTESTI*HLT 
JMS RTESTN 

11*0 

ISZ RTESTI 
JMP I RTESTI 


amax* 0; 0;0 
/ 

PAGE 

/ 

/ROUTINES 

RSAVAA# HLT /CALCULE : HOLU= (I PTAIJ) ;ii KIA1J ) = 
ClF+F'PPFLD 

JMS 17/ (I PTAR) = HOt-iJ 

FGET I PTAIJ 

PPUT HOLD 

FGET 1 PTAR 

FMPY FM1 

FPUT I PTAIJ 

FGET HOLD 

FPUT I PTAR 

FEXT 

JMP I RSAVAA 

/ 

/CALCUL DE HOLD= (I PTAIJ) 

/ (I PTAIJ) = C I PTAR) )(I PTAR) =-HULD 

RHOLAA# HLT 

cif+fppfld 

JMS I 7 
FGET I PTAIJ 
FPUT HOLD 
FGET I PTAR 
FPUT I PTAIJ 
FGET HOLD 
FMPY FM1 
FPUT I PTAR 
FEXT 

JMP I RHOLAA 


ABSVL# HLT 
CIF+FPPFLD 
JMS I 7 
FGET I DAMAX 
FSGE /FAC>=0 ? 


-<i PTAR) 


12 



FMPY KM 1 /NON 
FPUT FTEMPI 
FGET I PTAR 
FSGE 

FMPY FM1 
FSUB FTEMP1 
FEXT 
TAU 45 

SPA CLA /t ABS(MAX) < UU = FAbS (ACi,J)> ? 

JMP /tjyWft i ABSVU /non 

ClF+FPPFLD 
JMS I 7 

FGET 1 PTAR /GUI * PRENDRE A<I,J> COMM E MAX 

FPUT I UAMAX 

FEXT 

TAD 1 Dll 

DCA I PT1K /RELtVER LES INDICES 
TAD 1 DJ 
DCA I PTJK 
JMP 1 ABSVL 

/ 

/ 

/ 

DETAMA>HLT 
ClF+FPPFLD 
JMS 1 7 
FGET FUN 
FDIV I DAMAX 
FPUT I PTAR / l./AMAX 
FGET 1 DDET 
FMPY I UAMAX 
FPUT i DDET 
FEXT 

JMP 1 DETAMA 

/ 

/ 

APAMA* HLT 
JMS I CAL IK 
TAD PTAIJ 
DCA PTAR 
JMS I CALKJ 
ClF+FPPFLD 
JMS 1 7 
FGET 1 PTAIJ 
FMPY I PTAR 
FEXT 

JMS I CAL1J 
ClF+FPPFLD 
JMS 1 7 
FADD 1 PTAR 

FPUT 1 PTAR /A<1 * J)=A(1# J > + A(1>K)*A(K»J) 

FEXT 

JMP I APAMA 
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/ 

HOLD* 0*0*0 
FM1*II6000*0 
FUN*1*2000*0 
FTEMP1 * 0* 0* 0 
CALIK*RCALIK 
CALiJ*RCALIJ 
CALKJ* RCALKJ 
DAMAX* ANAX 
DDET* det 
DU * 11 
DJ* J 
/ 

DEB1K=. 

/ 

PAUSE 
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